#########################################################################
# Time-stamp: <05/09/09 11:06:23 ostolop>
#########################################################################
package EP::EPC::Core::Chromosome_Localization;

use strict;

use EP::Config;
use EP::Common;
use EP::Common::General;
use EP::Common::Error qw (:try);

use EP::Rclient;

use EP::EPC::Core;
use EP::EPC::Core::Component;
use EP::EPC::Core::Chromosome;

use File::Temp;
use File::Basename;

use Math::Complex;

#use SVG;

#Misha's SVG lib
use XML::LibXML::SVG;

our @ISA = qw ( EP::EPC::Core::Component );

##########################################################
###   GLOBAL VARIABALES                              #####
##########################################################
#Due to how it was written the chromosome localization requires
#a number of global variables. Will talk to Misha about the
#implications of this.

my %bands;
my %bandCount;
my %bandProb;
my %composites;
my %featureLocations;
my $compCount = 0;
my $setsize;
my $sample;
my %bandNorms;
my @cache;
my $pi      = 3.14159265358979323846264;
my $root2pi = sqrt( 2 * $pi );
my @chromosomes;
my %positions;

#The hash table which stores the features as a hash of arrays. Added to allow popups with the features inside;
my %featurebands;
my $timeC = 0;

#Mods to build  hashes from properties file
my %ideograms;
my %colours;

my %fullMembers;

#my %selectedMembers; This feature is currently being suppressed
my %bandOutput;
my $mostClones = 0;
my $species;    #This one should definitely be deglobalized
my $baseRefDir = $EP::Config::h->{CHROMOSOME_LOCALIZATION_IDEOGRAMS} . "/";
################################################################
#----------------------

sub process_component {
  my $self = shift;

  my ( $clg_table_url, $clg_image_url ) = $self->run_clg();

  my $process_result = {
    clg_table => {
      type  => "xml/table",
      value => $clg_table_url,
      label => "Chromosome Localization Probabilities"
    },
    clg_visualization => {
      type  => "image/png",
      value => $clg_image_url,
      label => "Chromosome Localization Plot"
    },
    clg_visualization_width  => 1500,
    clg_visualization_height => 700,

    setsize => {
      type  => "number",
      value => $setsize,
      label => "Set size"
    },
    samplesize => {
      type  => "number",
      value => $sample,
      label => "Sample size"
    },
    src_dataset_id => $self->{src_dataset}->id,
    dst_dataset_id => $self->{src_dataset}->id
  };

  $self->{CD}->{process_result} = $process_result;

  return $self->record_analysis_history();
}

#######################################################################
# Sub: run_clg
#   Runs the chromosome localization calculations
# Returns:
#   URL to an XML file containing the localization calculations table,
#   Name of the SVG file containing the localization graph
######################################################################
sub run_clg {
  my $self = shift;

  throw EP::Common::Error( -value => AFFY_TRANSFORM_RAW_DATA )
    if $self->{src_dataset}->type eq "affy";

  my $dataset_folder_path =
    $EP::Config::h->{EPNG_USRDATA_PATH} . "/" . $self->{src_dataset}->folder->location;

  ####
  # Read in parameters from cgi call

  my $query   = $self->{query};
  my $probMax = $query->param('probMax');
  my $percMin = $query->param('percMin');
  my $regMin  = $query->param('regMin');
  $species = $query->param('species');

  my $cloneLibrary = $query->param('cloneLibrary');
  my $filter       = $query->param('filter');
  if ( $filter eq "on" ) {
    $filter = 1;
  }
  my $affy_chip_id = $query->param('affy_chip_id');
  my $numericGraph = $query->param('numericGraph');
  if ( $numericGraph eq "on" ) {
    $numericGraph = 1;
  }
  my $numericMax = $query->param('numericMax');
  my $allDetails = $query->param('allD');
  if ( $allDetails eq "on" ) {
    $allDetails = 1;
  }
  my $percentage = $query->param('percent');
  if ( $percentage eq "on" ) {
    $percentage = 1;
  }
  my $maxper = $query->param('maxper');

  #Added For Bonferroni Correction
  my $bonferroni = $query->param('bonferroni');
  if ( $bonferroni eq "on" ) {
    $bonferroni = 1;
  } else {
    $bonferroni = 0;
  }
  ###
  # Initialize hash tables from properties file
  _load_properties();

  #####################################################################################################
  # Need to define the source dataset before the reference library so that the root dataset can be used
  # as a reference file.
  #####################################################################################################

  my $src_dataset  = $self->{src_dataset};
  my $root_dataset = $self->{session}->current_root_dataset;

  # check species karyotype matches the loaded dataset for datasets with known species
  # This has been generating a problem with the selections have to make sure the
  # dataset is taken from the root dataset as opposed to the src dataset. Is it likely
  # that the species will change through the session.

  if ( $root_dataset->species ne "unknown" ) {
    die "Karyotype species (hsapiens) and dataset species ("
      . $root_dataset->species
      . ") do not match."
      if $species =~ /hsapiens/ && $root_dataset->species ne "hsapiens";

    die "Karyotype species (mmusculus) and dataset species ("
      . $root_dataset->species
      . ") do not match."
      if $species =~ /mmusculus/ && $root_dataset->species ne "mmusculus";

    die "Karyotype species (dmelanogaster) and dataset species ("
      . $root_dataset->species
      . ") do not match."
      if $species =~ /dmelanogaster/ && $root_dataset->species ne "dmelanogaster";

    die "Karyotype species (rnorvegicus) and dataset species (" . $root_dataset->species . ") do not match."
      if $species =~ /rnorvegicus/ && $root_dataset->species ne "rnorvegicus";
  }

  #####################################################################################################
  # We now have a change in the method of loading the data for the reference file
  # Data can be taken from a number of sources namely custom reference, affymetrix id, arrayexpress and
  # the default is the root dataset for the sample being analyzed.
  # The processing for each of these sources is somewhat different and so they are dealt with in individual
  # subroutines.
  #####################################################################################################

  my $fullfile = $query->upload("reference_localization");
  if ($fullfile) {
    # A composite reference file is present this should be used as the population set
    countCompositesFromRefFH($fullfile);
  } elsif ( ($affy_chip_id) && $affy_chip_id ne 'none' ) {
    die "This functionality is temporarily disabled";
    # biomart retrieval
    EP::Rclient->R(<<__RCODE__);
app_root = "$EP::Config::h->{ROOT}";
library(biomaRt);
library(affy);
    
mart = martConnect();

cdf = cdfFromBioC( "$affy_chip_id", lib = "$EP::Config::h->{RLIB}/" );
probeIds = ls(cdf);

g = getGene( probeIds, array = "$affy_chip_id", mart = mart );
t = data.frame( feature_id = g\@id
              , chromosome_location = paste( g\@table\$chromosome, g\@table\$band, sep = "" )
              , composite_feature_id = g\@table\$symbol
              );
write.table( t, file="$EP::Config::h->{TMP}\/$affy_chip_id.chrocoloc", quote = FALSE, row.names = FALSE, sep = "\t" );
__RCODE__

    my $affy_chip_id_fh;
    open $affy_chip_id_fh, "$EP::Config::h->{TMP}/$affy_chip_id.chrocoloc"
      or die "can't open $EP::Config::h->{TMP}/$affy_chip_id.chrocoloc";
    countCompositesFromRefFH($affy_chip_id_fh);
    close $affy_chip_id_fh;
  } else {
    # Use the root dataset as the default. This should be representative of the original population tested
    my @root_ref_locs = map {
      [
        $root_dataset->data->rows->{$_}{full_annot}
          ->Values( $root_dataset->data->rows->{$_}{full_annot}->Indices('feature_id') ),
        $root_dataset->data->rows->{$_}{full_annot}
          ->Values( $root_dataset->data->rows->{$_}{full_annot}->Indices('chromosome_location') ),
        $root_dataset->data->rows->{$_}{full_annot}
          ->Values( $root_dataset->data->rows->{$_}{full_annot}->Indices('composite_feature_id') )
      ]
      } sort {
      $a <=> $b
      } keys %{ $root_dataset->data->rows };
    warn "GOING WITH INTERNAL ANNOTATION ", scalar @root_ref_locs, "\n";
    countCompositesFromRootRefArr( \@root_ref_locs );
  }

  # Reference lib should now be loaded in memory
  ####################################################################################################
  ### GETTING ANNOTATION COLUMNS ###
  
  # will try to use index of the 'feature_id' column if it's there, otherwise just use the first column (index=0)
  my $feature_id_idx = $self->{src_dataset}->data->rows->{0}{full_annot}->Indices('feature_id');
  $feature_id_idx = 0 if not defined $feature_id_idx;
  warn "Using $feature_id_idx for annotations";
  
  my @experimental_locs = map {
    [ $self->{src_dataset}->data->rows->{$_}{full_annot}
        ->Values( $feature_id_idx ) ]
    } sort {
    $a <=> $b
    } keys %{ $self->{src_dataset}->data->{row} };

  countCompositesPerBandExp( \@experimental_locs );
  ###########################################################################
  ###At this point all the initialization of the data sets should be done ###
  ###########################################################################

  ### DO YOUR CALCULATIONS ###

  #bandCount is now calculated
  #All band scores are put in a hash bandOutput as an array of results
  my @bks = keys %bandCount;
  calcBandStats( $filter, $probMax, $regMin, $percMin, $bonferroni );

#------------------
# place outputs into the folder in $dataset_folder_path
# should be two of these: the table and the svg image.
# the URL to the table is formed by taking $output_www_base (below) and the name of the file you write the data to. see below for example.

  my $output_www_base =
    "$EP::Config::h->{EPNG_USRDATA_URL}/" . $self->{src_dataset}->folder->location;

  my $clg_table_file = new File::Temp(
    DIR      => $dataset_folder_path,
    TEMPLATE => "clg_table_XXXXXX",
    SUFFIX   => ".xml",
    UNLINK   => 0
  );
  my $clg_table_filename = basename($clg_table_file);
  my $clg_table_url      = $output_www_base . "/" . $clg_table_filename;

  my $clg_svg_file = new File::Temp(
    DIR      => $dataset_folder_path,
    TEMPLATE => "clg_svg_XXXXXX",
    SUFFIX   => ".svg",
    UNLINK   => 0
  );
  my $clg_svg_filename   = basename($clg_svg_file);
  my $clg_svg_url        = $output_www_base . "/" . $clg_svg_filename;

  my $xml = qq{<?xml version='1.0'?>};
  $xml .= "<table>";

  $xml .= "<head>";
  $xml .= "<data>$_</data>"
    foreach ( "Chromosome Band", "Regulated per Band", "Present per Band", "Reg/Present",
    "P-value" );
  $xml .= "</head>";

  my @bandoutkeys = keys %bandOutput;
  my $sorted      = sortExpKeys( \@bandoutkeys );
  foreach my $boKey (@$sorted) {
    $xml .= "<row>";
    my $row = $bandOutput{$boKey};
    $xml .= "<data>$_</data>" foreach @$row;
    $xml .= "</row>";
  }

  $xml .= "</table>";

  print $clg_table_file $xml;
  close $clg_table_file;

  # CREATE THE SVG HERE, WRITE TO a file accessible through URL $output_www_base . $clg_svg_filename
  #Additions for SVG
  my %colouredBands;
  if ( $numericGraph == 1 ) {
    if ( $numericMax > 0 ) {
      $mostClones = $numericMax;
    }
    $probMax = $mostClones;
    calcSFshashNumeric( \%colouredBands, $mostClones );
  } elsif ( $percentage == 1 ) {
    calcSFshashPercentage( \%colouredBands, $maxper );
    $probMax = $maxper;
  } else {

    calcSFshash( \%colouredBands, $probMax, $regMin, $percMin );
  }

  my $svg_xml = _SVGccbands( \%colouredBands, $probMax, $numericGraph, $percentage );
  print $clg_svg_file $svg_xml;
  close $clg_svg_file;

  # convert SVG to PNG
  warn "Rasterizing $clg_svg_file ...";

  my $clg_image_file = $clg_svg_file;
  $clg_image_file =~ s/\.svg/\.png/;

  if ( system_call("$EP::Config::h->{CONVERT_CMD} $clg_svg_file $clg_image_file") == 0 ) {
    die "Failed to convert $clg_svg_file to $clg_image_file.";
  }

  my $clg_image_url = $clg_svg_url;
  $clg_image_url =~ s/\.svg/\.png/;

  return ( $clg_table_url, $clg_image_url );
}
#############################
# Additions by Jonathon Blake
#############################
###################################################
# SUBROUTINES
###################################################

####################################################################################################
#Title: countCompositesPerBandExp
#Description: Take an array of ids and checks against the chromosome banding via composite ids and
#             featureLocations populated earlier
#Author: Jon
#Date: 22.07.2005
####################################################################################################
sub countCompositesPerBandExp {
  my ($ids) = @_;

  #local composite control
  my %comps;

  # print STDERR "THERE ARE ", scalar @$ids, " IDS\n";
  foreach my $id (@$ids) {

    if ( exists $composites{ $id->[0] } ) {

      #feature has a chromosome location only those that do are loaded from reference file
      #Has the composite been seen before
      if ( !( exists $comps{ $composites{ $id->[0] } } ) ) {
        $bandCount{ $featureLocations{ $composites{ $id->[0] } } }++;
        $sample++;
        $comps{ $composites{ $id->[0] } } = 1;

        #Additions to allow the production of the popup windows with features
        #Check to see if there is already an array for the feature in the bandfeatures hash
        if ( exists $featurebands{ $featureLocations{ $composites{ $id->[0] } } } ) {
          push( @{ $featurebands{ $featureLocations{ $composites{ $id->[0] } } } }, $id->[0] );
        } else {

          #Array does not exist create one
          my @featureArr;
          push( @featureArr, $id->[0] );
          $featurebands{ $featureLocations{ $composites{ $id->[0] } } } = \@featureArr;
        }
      }
    }
  }

  undef %comps;
}
####################################################################################################
#-------------------------
sub fac {
  my $n = shift;

  # testing
  return $cache[$n] if defined $cache[$n];

  # testing
  return undef if $n < 0;
  return 1     if $n == 0;
  return $cache[$n] = $n * fac( $n - 1 );
}

#------------------------
#clusterProbability
# calculates the probability of a cluster of genes according to the method of Chang et al

sub clusterProbability {
  my ( $x, $y, $z, $win ) = @_;    # Variables named after Chang et al
                                   #A safety net build in to stop log of 0 errors
       #Taking that in most circumstances the probability of all the set being up upregulated is small.
       # return zero if the cluster equals the window size
  my $sfwCg = 0;
  if ( $z == $win ) {

    #return 0;
    #need a special value for wCg taking into account that factorial 0 = 1
    facDurchfac( $win, $z );

  } else {
    $sfwCg = xCy( $win, $z );
  }
  my $sfnCr = xCy( $x, $y );
  my $sfCombCor = xCy( ( $x - $win ), ( $y - $z ) );
  my $calcApprox = $sfwCg + ( $sfCombCor - $sfnCr );
  my $approxPro = exp($calcApprox);
  return $approxPro;
}

#------------------------------------------------------------------
sub facDurchfac {
  my ( $n, $x ) = @_;

  #returns n!/x!
  #Use the Schwager simplification exp((n+1/2)*ln(n) -((x+1/2)*ln(x))+x -n)
  return ( ( $n + 0.5 ) * log($n) - ( ( $x + 0.5 ) * log($x) ) + $x - $n );
}

#------------------------------------------
#-------------------------------------------
# calcBandStats
# Calculates the hypogeometric distribution for the bands in comparision to the genome.
sub calcBandStats {
  my ( $filter, $maxprob, $regmin, $percmin, $bonferroni ) = @_;
  my @expKeys = keys %bandCount;

#%bandCount is built from the sample data and should represent all the windows ie tests used in the sampling
#Therefore I am using this as the basis of the Bonferroni calculation
  my $bonCor = scalar @expKeys;
  print STDERR "The bonferroni = $bonCor\n";

  #sort keys
  my $sortedKeys = sortExpKeys( \@expKeys );
  if ( $filter == 1 ) {
    foreach my $expkey (@$sortedKeys) {

      #Whole genome exp probability. This is the correct one
      my $prob = clusterProbability( $setsize, $sample, $bandCount{$expkey}, $bands{$expkey} );
      ###
      #Bonferroni Correction for Multiple testing
      if ( $bonferroni == 1 ) {
        $prob *= $bonCor;
        if ( $prob > 1 ) {
          $prob = 1;
        }
      }
      ###
      if ( $bandCount{$expkey} > $mostClones ) {
        $mostClones = $bandCount{$expkey};
      }

      #At this point we could see if we have a hotspot, a coldspot or nothing
      my $bandNorm = $bandCount{$expkey} / $bands{$expkey};

      #Apply filter to output
      $bandProb{$expkey}  = $prob;
      $bandNorms{$expkey} = $bandNorm;
      if ( ( $bandCount{$expkey} > $regmin && $bandNorm >= $percmin / 100 ) && $prob <= $maxprob ) {

        $bandProb{$expkey}  = $prob;
        $bandNorms{$expkey} = $bandNorm;
        my @tabrow = (
          $expkey, $bandCount{$expkey}, $bands{$expkey},
          sprintf( "%.4f", $bandNorm ),
          sprintf( "%.4f", $prob )
        );
        $bandOutput{$expkey} = \@tabrow;
      }
    }

  } else {
    foreach my $expkey (@$sortedKeys) {
      my $prob = clusterProbability( $setsize, $sample, $bandCount{$expkey}, $bands{$expkey} );
      ###
      #Bonferroni Correction for Multiple testing
      if ( $bonferroni == 1 ) {
        $prob *= $bonCor;
        if ( $prob > 1 ) {
          $prob = 1;
        }
      }
      ###
      #At this point we could see if we have a hotspot, a coldspot or nothing
      my $bandNorm = $bandCount{$expkey} / $bands{$expkey};
      $bandProb{$expkey}  = $prob;
      $bandNorms{$expkey} = $bandNorm;
      my @tabrow = (
        $expkey, $bandCount{$expkey}, $bands{$expkey},
        sprintf( "%.4f", $bandNorm ),
        sprintf( "%.4g", $prob )
      );
      $bandOutput{$expkey} = \@tabrow;
    }
  }

}
############################################################################
############################################################################
sub xCy {
  my ( $x, $y ) = @_;
  my $lnFy = 0;

  #First build the logs
  #exp of x! for division
  my $lnFx = ( log($x) * ( $x + 0.5 ) ) - $x;

  #Add a correction for when $x = $y
  my $lnFxy;
  if ( $x == $y ) {
    $lnFxy = 0;
  } else {
    $lnFxy = ( log( $x - $y ) * ( ( $x - $y ) + 0.5 ) - ( $x - $y ) );
  }

  #Taking the following as attempt to increase accurracy
  if ( $y > 20 ) {
    $lnFy = log($root2pi) + ( log($y) * ( $y + 0.5 ) ) - $y;
  } else {
    $lnFy = log( fac($y) );
  }
  my $approx = $lnFx - ( $lnFxy + $lnFy );
  return ($approx);    # - (log($x-$y)*(($x-$y)+0.5)+(log($y)*($y+0.5)));
}

#------------------------------------
sub sort_orig_out_bandRef {
  my ($cband) = @_;

  #Sort out band.
  #Is it a band range
  my $chrom = "";
  $cband =~ s/\s*//;
  if ( $cband =~ /(.*)-(.*)/ ) {
    $cband = $1;
    my $band2 = $2;
    if ( $cband =~ /q11/ ) {
      $cband =~ /(\S*)q/;
      $cband = $1 . $band2;
    }
  }

  #If it is a sub band that is not found

  if ( $cband =~ /^(\S*?)[p|q]/ ) {
    $chrom = $1;
  }
  my @retarr = ( $cband, $chrom );
  return @retarr;

}

#------------------------
#Bits added from Chromosome Test handle
#_initChroms :- a new parser for the chromosome format makes a lot of chromosomes.
sub _initChroms {

  #Open ideogram file
  my $refile = $baseRefDir . $ideograms{$species};
  open( CIN, $refile ) or die "Can't open cbands $refile\n";
  my $chromcounter = 0;
  my $currentChrom = 0;
  my $qarm;
  my $parm;
  my $bandmode;
  my $acro     = "no";
  my $chromref = \@chromosomes;

  while (<CIN>) {
    my $linedata = $_;
    $linedata =~ s/[\r|\b]//ig;
    if ( $linedata =~ /CH:\"(\S*)\";P:\"(\d*\.\d*)\";Q:\"(\d*\.\d*)\";/ ) {
      my $chrom = new Chromosome( $1, $2, $3 );
      $chromosomes[$chromcounter] = $chrom;
      $chromcounter++;
      $currentChrom = $chromcounter - 1;
      $chromosomes[$currentChrom]->{acro} = "no";
    } elsif ( $linedata =~ /qbands;/ ) {
      $bandmode = "q";
    } elsif ( $linedata =~ /pbands;/ ) {
      $bandmode = "p";
    } elsif (
      $linedata =~ /name:\"(\S*)\";.*colour:(\S*);.*stop:\"(\d*\.\d*)\";.*start:\"(0\.\d*)\";/ )
    {
      $chromosomes[$currentChrom]->add_Band( $bandmode, $1, $2, $4, $3 );
    }

    #Need an extra line to make sure we get the last band (No decimal) in
    elsif ( $linedata =~ /name:\"(\S*)\";.*colour:(\S*);.*stop:\"(\d*)\";.*start:\"(0\.\d*)\";/ ) {
      $chromosomes[$currentChrom]->add_Band( $bandmode, $1, $2, $4, $3 );
    }

    #Need an extra line to make sure we get the first band(No decimal) in
    elsif ( $linedata =~ /name:\"(\S*)\";.*colour:(\S*);.*stop:\"(\d*\.\d*)\";.*start:\"(0)\";/ ) {
      $chromosomes[$currentChrom]->add_Band( $bandmode, $1, $2, $4, $3 );
    } elsif ( $linedata =~ /ACROCENTRIC/ ) {
      $chromosomes[$currentChrom]->{acro} = "yes";
    }
  }
  close(CIN);

}
##################################################################################################

sub _init_positions {
  my ( $h, $w, $n ) = @_;

  #Top margin is 10% of total height
  my $c = 0;

  #Drosophila has the 4th chromosome as the largest
  if ( $species eq "drosophila" ) {
    $c = 3;
  }
  my $tm = $h * 0.1;

  #my %positions = $pos;
  # The drawable face of the paper is 60% of total height
  $positions{height} = $h * 0.6;

  #The centromere position is calculated using the largest chromosome
  my $parmPecentage =
    $chromosomes[$c]->{parm} / ( $chromosomes[$c]->{parm} + $chromosomes[$c]->{qarm} );
  my $qarmPecentage =
    $chromosomes[$c]->{qarm} / ( $chromosomes[$c]->{parm} + $chromosomes[$c]->{qarm} );
  $positions{parm}  = $parmPecentage;
  $positions{qarm}  = $qarmPecentage;
  $positions{ycent} = ( $positions{height} * $parmPecentage ) + $tm;

  #Width of the chromosome = 90% of the width of paper / chrom number *2
  $positions{width} = ( $w * 0.9 ) / ( $n * 2 );

  #Puts a maximum with on the chromosomes
  if ( $positions{width} > 50 ) {
    $positions{width} = 50;
  }

  #First margin is 5% of total width;
  $positions{xcent} = $w * 0.05;
}
################################################################################################
sub _SVGccbands {

  #Receive a reference to the coloured band hash on calling.
  my ( $colouredBands, $maxprob, $nummericGraph, $percentage ) = @_;
  _initChroms();

  #As this is for a server and no configuration file is currently present.
  #HARDCODE of height and width follows.
  my $height = 700;
  my $width  = 1500;
  _init_positions( $height, $width, scalar @chromosomes );

  #--------------------
  #Misha wants to use his SVG module, so changes added here.
  #EPC/Visualization has been used as a basis.
  # my $svg = new SVG->new(width=>$width, height=>$height);
  my $svg_params = {
    'xmlns:a3'                => "http://ns.adobe.com/AdobeSVGViewerExtensions/3.0/",
    'xmlns:xlink'             => "http://www.w3.org/1999/xlink",
    'xmlns:gui'               => "http://www.kevlindev.com/gui",
    'a3:scriptImplementation' => "Adobe",
    preserveAspectRatio       => "xMinYMin",
    onzoom                    => "on_zoom()",

    #  onload                    => "init(evt)",
    width            => $width,
    height           => $height,
    viewBox          => "0 0 $width $height",
    'text-rendering' => 'optimizeSpeed',

    #                    'shape-rendering'          => 'optimizeSpeed',
    'image-rendering' => 'optimizeSpeed',
    'color-rendering' => 'optimizeSpeed'
  };

  my $svg = XML::LibXML::SVG->createSVGDocument(%$svg_params);

  #--------------------
  my $c = 0;
  if ( $species eq "drosophila" ) {
    $c = 3;
  }
  my $parmPc = $chromosomes[0]->{parm} / ( $chromosomes[0]->{parm} + $chromosomes[$c]->{qarm} );
  my $qarmPc = $chromosomes[0]->{qarm} / ( $chromosomes[0]->{parm} + $chromosomes[$c]->{qarm} );
  $positions{qarm} = $qarmPc * ( $chromosomes[0]->{qarm} / $chromosomes[$c]->{qarm} );

  $chromosomes[0]->draw_self( $svg, $colouredBands, \%positions, \%featurebands );

  #For each of the following chromosomes
  for ( my $x = 1 ; $x < scalar @chromosomes ; $x++ ) {
    $positions{xcent} += 2 * ( $positions{width} );

    #Change for Drosophila
    if ( $chromosomes[$x]->{parm} > 0 ) {
      $positions{parm} = $parmPc * ( $chromosomes[$x]->{parm} / $chromosomes[$c]->{parm} );
    }
    $positions{qarm} = $qarmPc * ( $chromosomes[$x]->{qarm} / $chromosomes[$c]->{qarm} );

    $chromosomes[$x]->draw_self( $svg, $colouredBands, \%positions, \%featurebands );
  }

  if ( $nummericGraph == 1 ) {
    my %cStyle = ( 'stroke' => 'black', 'fill' => $colours{RED}, 'stroke-width' => '0' );
    my $legend = "> " . ($maxprob);
    $svg->rect(
      x      => $width * 0.1,
      y      => $height * 0.9,
      width  => 20,
      height => 10,
      id     => 'col5',
      %cStyle
    );
    $svg->text(
      id    => 'col5label',
      x     => ( ( $width * 0.1 ) + 25 ),
      y     => ( $height * 0.9 ) + 10,
      style => { 'font' => 'Arial', 'font-size' => 10, 'fill' => 'black' }
    )->cdata($legend);
    $legend = "<= " . ($maxprob);
    %cStyle = ( 'stroke' => 'black', 'fill' => $colours{ORANGE}, 'stroke-width' => '0' );
    $svg->rect(
      x      => $width * 0.1,
      y      => ( $height * 0.9 ) + 10,
      width  => 20,
      height => 10,
      id     => 'col4',
      %cStyle
    );
    $svg->text(
      id    => 'col4label',
      x     => ( ( $width * 0.1 ) + 25 ),
      y     => ( $height * 0.9 ) + 20,
      style => { 'font' => 'Arial', 'font-size' => 10, 'fill' => 'black' }
    )->cdata($legend);
    $legend = "<= " . ( 3 * ( $maxprob / 4 ) );
    %cStyle = ( 'stroke' => 'black', 'fill' => $colours{YELLOW}, 'stroke-width' => '0' );
    $svg->rect(
      x      => $width * 0.1,
      y      => ( $height * 0.9 ) + 20,
      width  => 20,
      height => 10,
      id     => 'col3',
      %cStyle
    );
    $svg->text(
      id    => 'col3label',
      x     => ( ( $width * 0.1 ) + 25 ),
      y     => ( $height * 0.9 ) + 30,
      style => { 'font' => 'Arial', 'font-size' => 10, 'fill' => 'black' }
    )->cdata($legend);
    $legend = "<= " . ( 2 * ( $maxprob / 4 ) );
    %cStyle = ( 'stroke' => 'black', 'fill' => $colours{GREEN}, 'stroke-width' => '0' );
    $svg->rect(
      x      => $width * 0.1,
      y      => ( $height * 0.9 ) + 30,
      width  => 20,
      height => 10,
      id     => 'col2',
      %cStyle
    );
    $svg->text(
      id    => 'col2label',
      x     => ( ( $width * 0.1 ) + 25 ),
      y     => ( $height * 0.9 ) + 40,
      style => { 'font' => 'Arial', 'font-size' => 10, 'fill' => 'black' }
    )->cdata($legend);
    $legend = "<= " . ( $maxprob / 4 );
    %cStyle = ( 'stroke' => 'black', 'fill' => $colours{LIGHTGREEN}, 'stroke-width' => '0' );
    $svg->rect(
      x      => $width * 0.1,
      y      => ( $height * 0.9 ) + 40,
      width  => 20,
      height => 10,
      id     => 'col1',
      %cStyle
    );
    $svg->text(
      id    => 'col1label',
      x     => ( ( $width * 0.1 ) + 25 ),
      y     => ( $height * 0.9 ) + 50,
      style => { 'font' => 'Arial', 'font-size' => 10, 'fill' => 'black' }
    )->cdata($legend);
  } elsif ( $percentage == 1 ) {
    my %cStyle = ( 'stroke' => 'black', 'fill' => $colours{RED}, 'stroke-width' => '0' );
    my $legend = "> " . $maxprob . "%";
    $svg->rect(
      x      => $width * 0.1,
      y      => $height * 0.9,
      width  => 20,
      height => 10,
      id     => 'col5',
      %cStyle
    );
    $svg->text(
      id          => 'col5label',
      x           => ( ( $width * 0.1 ) + 25 ),
      y           => ( $height * 0.9 ) + 10,
      'font-size' => 10,
      'fill'      => 'black'
    )->cdata($legend);
    $legend = "<= " . $maxprob . "%";
    %cStyle = ( 'stroke' => 'black', 'fill' => $colours{ORANGE}, 'stroke-width' => '0' );
    $svg->rect(
      x      => $width * 0.1,
      y      => ( $height * 0.9 ) + 10,
      width  => 20,
      height => 10,
      id     => 'col4',
      %cStyle
    );
    $svg->text(
      id    => 'col4label',
      x     => ( ( $width * 0.1 ) + 25 ),
      y     => ( $height * 0.9 ) + 20,
      style => { 'font' => 'Arial', 'font-size' => 10, 'fill' => 'black' }
    )->cdata($legend);
    $legend = "<= " . ( 3 * ( $maxprob / 4 ) ) . "%";
    %cStyle = ( 'stroke' => 'black', 'fill' => $colours{YELLOW}, 'stroke-width' => '0' );
    $svg->rect(
      x      => $width * 0.1,
      y      => ( $height * 0.9 ) + 20,
      width  => 20,
      height => 10,
      id     => 'col3',
      %cStyle
    );
    $svg->text(
      id    => 'col3label',
      x     => ( ( $width * 0.1 ) + 25 ),
      y     => ( $height * 0.9 ) + 30,
      style => { 'font' => 'Arial', 'font-size' => 10, 'fill' => 'black' }
    )->cdata($legend);
    $legend = "<= " . ( 2 * ( $maxprob / 4 ) ) . "%";
    %cStyle = ( 'stroke' => 'black', 'fill' => $colours{GREEN}, 'stroke-width' => '0' );
    $svg->rect(
      x      => $width * 0.1,
      y      => ( $height * 0.9 ) + 30,
      width  => 20,
      height => 10,
      id     => 'col2',
      %cStyle
    );
    $svg->text(
      id    => 'col2label',
      x     => ( ( $width * 0.1 ) + 25 ),
      y     => ( $height * 0.9 ) + 40,
      style => { 'font' => 'Arial', 'font-size' => 10, 'fill' => 'black' }
    )->cdata($legend);
    $legend = "<= " . ( $maxprob / 4 ) . "%";
    %cStyle = ( 'stroke' => 'black', 'fill' => $colours{LIGHTGREEN}, 'stroke-width' => '0' );
    $svg->rect(
      x      => $width * 0.1,
      y      => ( $height * 0.9 ) + 40,
      width  => 20,
      height => 10,
      id     => 'col1',
      %cStyle
    );
    $svg->text(
      id    => 'col1label',
      x     => ( ( $width * 0.1 ) + 25 ),
      y     => ( $height * 0.9 ) + 50,
      style => { 'font' => 'Arial', 'font-size' => 10, 'fill' => 'black' }
    )->cdata($legend);
  } else {
    my %cStyle = ( 'stroke' => 'black', 'fill' => $colours{RED}, 'stroke-width' => '0' );
    my $legend = "<= " . ( $maxprob / 5 );
    $svg->rect(
      x      => $width * 0.1,
      y      => $height * 0.9,
      width  => 20,
      height => 10,
      id     => 'col5',
      %cStyle
    );
    $svg->text(
      id    => 'col5label',
      x     => ( ( $width * 0.1 ) + 25 ),
      y     => ( $height * 0.9 ) + 10,
      style => { 'font' => 'Arial', 'font-size' => 10, 'fill' => 'black' }
    )->cdata($legend);
    $legend = "<= " . ( 2 * ( $maxprob / 5 ) );
    %cStyle = ( 'stroke' => 'black', 'fill' => $colours{ORANGE}, 'stroke-width' => '0' );
    $svg->rect(
      x      => $width * 0.1,
      y      => ( $height * 0.9 ) + 10,
      width  => 20,
      height => 10,
      id     => 'col4',
      %cStyle
    );
    $svg->text(
      id    => 'col4label',
      x     => ( ( $width * 0.1 ) + 25 ),
      y     => ( $height * 0.9 ) + 20,
      style => { 'font' => 'Arial', 'font-size' => 10, 'fill' => 'black' }
    )->cdata($legend);
    $legend = "<= " . ( 3 * ( $maxprob / 5 ) );
    %cStyle = ( 'stroke' => 'black', 'fill' => $colours{YELLOW}, 'stroke-width' => '0' );
    $svg->rect(
      x      => $width * 0.1,
      y      => ( $height * 0.9 ) + 20,
      width  => 20,
      height => 10,
      id     => 'col3',
      %cStyle
    );
    $svg->text(
      id    => 'col3label',
      x     => ( ( $width * 0.1 ) + 25 ),
      y     => ( $height * 0.9 ) + 30,
      style => { 'font' => 'Arial', 'font-size' => 10, 'fill' => 'black' }
    )->cdata($legend);
    $legend = "<= " . ( 4 * ( $maxprob / 5 ) );
    %cStyle = ( 'stroke' => 'black', 'fill' => $colours{GREEN}, 'stroke-width' => '0' );
    $svg->rect(
      x      => $width * 0.1,
      y      => ( $height * 0.9 ) + 30,
      width  => 20,
      height => 10,
      id     => 'col2',
      %cStyle
    );
    $svg->text(
      id    => 'col2label',
      x     => ( ( $width * 0.1 ) + 25 ),
      y     => ( $height * 0.9 ) + 40,
      style => { 'font' => 'Arial', 'font-size' => 10, 'fill' => 'black' }
    )->cdata($legend);
    $legend = "<= " . $maxprob;
    %cStyle = ( 'stroke' => 'black', 'fill' => $colours{LIGHTGREEN}, 'stroke-width' => '0' );
    $svg->rect(
      x      => $width * 0.1,
      y      => ( $height * 0.9 ) + 40,
      width  => 20,
      height => 10,
      id     => 'col1',
      %cStyle
    );
    $svg->text(
      id    => 'col1label',
      x     => ( ( $width * 0.1 ) + 25 ),
      y     => ( $height * 0.9 ) + 50,
      style => { 'font' => 'Arial', 'font-size' => 10, 'fill' => 'black' }
    )->cdata($legend);

  }

  $svg->group( id => "featureHook" );

  # my $xml = $svg->xmlify;
  my $xml = $svg->toString();
  return $xml;

}

#-------------------------
sub calcSFshash {
  my ( $sfbhash, $maxprob, $regmin, $percmin ) = @_;
  my @expKeys = keys %bandCount;

  foreach my $expkey (@expKeys) {
    my $bandVal   = $bandProb{$expkey};
    my $outString = "";
    if ( $bandCount{$expkey} > $regmin && $bandNorms{$expkey} >= $percmin / 100 ) {
      if ( $bandVal <= ( $maxprob / 5 ) ) {
        $sfbhash->{$expkey} = $colours{RED};    #'Col20plus';

      } elsif ( $bandVal <= 2 * ( $maxprob / 5 ) ) {
        $sfbhash->{$expkey} = $colours{ORANGE};    #'ColTo20';

      } elsif ( $bandVal <= 3 * ( $maxprob / 5 ) ) {
        $sfbhash->{$expkey} = $colours{YELLOW};    #'ColTo15';

      } elsif ( $bandVal <= 4 * ( $maxprob / 5 ) ) {
        $sfbhash->{$expkey} = $colours{GREEN};     #'ColTo10';

      } elsif ( $bandVal <= $maxprob ) {
        $sfbhash->{$expkey} = $colours{LIGHTGREEN};    #'ColTo5';

      }
    }
  }
}

#-------------------------
sub calcSFshashNumeric {
  my ( $sfbhash, $most ) = @_;
  my @expKeys = keys %bandCount;

  foreach my $expkey (@expKeys) {
    my $bandVal   = $bandCount{$expkey};
    my $outString = "";
    if ( $bandVal <= ( $most / 4 ) ) {
      $sfbhash->{$expkey} = $colours{LIGHTGREEN};    #'Col20plus';

    } elsif ( $bandVal <= 2 * ( $most / 4 ) ) {
      $sfbhash->{$expkey} = $colours{GREEN};         #'ColTo20';

    } elsif ( $bandVal <= 3 * ( $most / 4 ) ) {
      $sfbhash->{$expkey} = $colours{YELLOW};        #'ColTo15';

    } elsif ( $bandVal <= 4 * ( $most / 4 ) ) {
      $sfbhash->{$expkey} = $colours{ORANGE};        #'ColTo10';

    } elsif ( $bandVal > $most ) {
      $sfbhash->{$expkey} = $colours{RED};           #'ColTo5';
    }
  }
}

#--------------------------
sub calcSFshashPercentage {
  my ( $sfbhash, $most ) = @_;
  my @expKeys = keys %bandCount;
  foreach my $expkey (@expKeys) {
    my $bandVal   = 100 * ( $bandCount{$expkey} / $bands{$expkey} );
    my $outString = "";

    #print "regmin = $regmin, percmin = $percmin\n";
    if ( $bandVal <= ( $most / 4 ) ) {
      $sfbhash->{$expkey} = $colours{LIGHTGREEN};    #'Col20plus';
    } elsif ( $bandVal <= ( 2 * ( $most / 4 ) ) ) {
      $sfbhash->{$expkey} = $colours{GREEN};         #'ColTo20';
    } elsif ( $bandVal <= ( 3 * ( $most / 4 ) ) ) {
      $sfbhash->{$expkey} = $colours{YELLOW};        #'ColTo15';
    } elsif ( $bandVal <= $most ) {
      $sfbhash->{$expkey} = $colours{ORANGE};        #'ColTo10';
    } elsif ( $bandVal > $most ) {
      $sfbhash->{$expkey} = $colours{RED};           #'ColTo5';
    }
  }
}

#--------------------------
sub sortExpKeys {
  my ($expKeys) = @_;

  #Arrays required
  my @singleFigureChroms;
  my @XYChroms;
  my @doubleFigureChroms;

  #Go through each of the keys and assign to the correct sorting array
  foreach my $exp (@$expKeys) {
    if ( $exp =~ /^\d[p|q]/ ) {
      push( @singleFigureChroms, $exp );
    } elsif ( $exp =~ /^\w[p|q]/ ) {
      push( @XYChroms, $exp );
    } else {
      push( @doubleFigureChroms, $exp );
    }
  }

  #chrom bands are assigned, now sort each array.
  my @sorted1digit = sort { $a cmp $b } @singleFigureChroms;
  my @sorted2digit = sort { $a cmp $b } @doubleFigureChroms;
  my @sortedXY     = sort { $a cmp $b } @XYChroms;

  #rebuild them back together.
  my @sortedAll;
  my $numberSorted = 0;
  foreach my $sort (@sorted1digit) {
    $sortedAll[$numberSorted] = $sort;
    $numberSorted++;
  }
  foreach my $sort (@sorted2digit) {
    $sortedAll[$numberSorted] = $sort;
    $numberSorted++;
  }
  foreach my $sort (@sortedXY) {
    $sortedAll[$numberSorted] = $sort;
    $numberSorted++;
  }
  return \@sortedAll;
}

#*******************************************************************
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
sub countCompositesPerBandRefFile {

  #In this setup the reffile is now a large string
  my ( $refString, $cmp ) = @_;
  my %comps;
  my @refs = split( "\n", $refString );
  foreach my $ref (@refs) {
    my @data = split( "\t", $ref );

    #
    $data[$cmp] =~ s/\r|\b//g;
    chomp( $data[$cmp] );
    my ( $band, $chrom ) = sort_orig_out_bandRef( $data[1] );

    if ( $band =~ /(^\S*)[q|p]/ ) {

      if ( $data[$cmp] ne "" ) {
        if ( exists $comps{ $data[$cmp] } ) {

          next;
        } else {
          $comps{ $data[$cmp] } = 1;

          #This should show where we are adding stuff to. Composites should only be counted once
          $bands{$band}++;

          $setsize++;
        }
      } else {
        $bands{$band}++;
        $setsize++;

      }

    }
  }

  #As cleaning up undef comps

  undef %comps;
}
######################################################################################################
#-------------------------
# Title: countCompositesPerRootDataRef
# Description: Produces the counts of features present for the reference population set from
# the EP style data array
# Date: 21.07.2005
# Author: Jon
######################################################################################################
sub countCompositesPerRootDataRef {
  my ( $featureArray, $cmp ) = @_;
  my %comps;    # Tempory hash to check composites

  foreach my $feature (@$featureArray) {
    my ( $band, $chrom ) = sort_orig_out_bandRef( $feature->[1] );
    if ( $band =~ /(^\S*)[q|p]/ ) {
      if ( $feature->[$cmp] ne "" ) {
        if ( exists $comps{ $feature->[$cmp] } ) {
          next;
        } else {
          $comps{ $feature->[$cmp] } = 1;

          #This should show where we are adding stuff to. Composites should only be counted once
          $bands{$band}++;
          $setsize++;
        }
      } else {
        $bands{$band}++;
        $setsize++;

      }
    }
  }    #End of foreach loop
  undef %comps;
}
#####################################################################################################
#Title: countCompositesFromRefFH
#Description: Takes a reference file handle and populates the band counts for the reference sets
#             and the annotation for the sample set. composites and featureLocations
#Author: Jon
#Date: 22.07.2005
#####################################################################################################
sub countCompositesFromRefFH {
  my ($fh) = @_;

  #Got the filehandle, get and possibly check the first line
  my $linedata = <$fh>;
  print STDERR "Line DATA is $linedata\n";
  if ( !( $linedata =~ /^feature_id/ ) ) { die "Poorly formed reference file\n"; }

  #Assume not that every thing is okay
  while (<$fh>) {
    my @data = split( "\t", $_ );
    chomp( $data[2] );    #The 3rd column is likely to be the last, but it won't die if it isn't
                          #Decide on the quality of chromosome assignment
                          #e.g is it there
    if ( ( !( $data[1] ) ) || $data[1] =~ /NA/i ) {

#Not wasting memory on features with no details
#Not risking that a composite with  data on one line is blocked by an earlier entry without annotation
      next;
    }
    if ( length( $data[1] ) < 4 ) {

#Chromosome annotation varies greatly across species, which is a pain and stops one doing a proper check
#All human, rat etc meaningful annotations will be longer than 3 characters e.g 1p21 so this is the simplest filter
#I can put into block chromosome or arm only annotations.
      next;
    }

    #Occassionally need to do some tidying up of the chromosome location
    my $cloc;
    if ( $data[1] =~ /(.*)-.*/ ) {
      $cloc = $1;
    } else {
      $cloc = $data[1];
    }
    $cloc =~ s/\s*//ig;    #have to trim these things
                           #Load the feature into the composite and featureLocations hashes
    if  (( !( $data[2] ) ) || $data[2] =~ /NA/ ) {

      #There is no composite here so the setsize gets incremented
      $setsize++;
      $composites{ $data[0] }       = $data[0];    # The cloneid is its own composite id
      $featureLocations{ $data[0] } =
        $cloc;            #This effectively using the composite id as they are the same in this case
      $bands{$cloc}++;    #increment the band/window count
    } else {

      #A composite is present
      $composites{ $data[0] } = $data[2];

      #Check to see if it has not been already counted
      if ( !( exists $featureLocations{ $data[2] } ) ) {

        #The composite has not yet been added, update and increment counts
        $featureLocations{ $data[2] } = $cloc;
        $setsize++;
        $bands{$cloc}++;
      }
    }
  }
}

#-------------------------
#####################################################################################################
#Title: countCompositesFromRefArr
#Description: Takes an array of data from the root data set and populates the band counts for the reference sets
#             and the annotation for the sample set. composites and featureLocations
#Author: Jon
#Date: 22.07.2005
#####################################################################################################
sub countCompositesFromRootRefArr {
  my ($rows) = @_;
  foreach my $data (@$rows) {

    #Decide on the quality of chromosome assignment
    #e.g is it there
   # if ( ( !( $data->[1] ) ) || $data->[1] eq "NA" ) {
     if  (( !( $data->[1] ) ) || $data->[1] =~ /NA/i ) {

#Not wasting memory on features with no details
#Not risking that a composite with  data on one line is blocked by an earlier entry without annotation
      next;
    }
    if ( length( $data->[1] ) < 4 ) {

#Chromosome annotation varies greatly across species, which is a pain and stops one doing a proper check
#All human, rat etc meaningful annotations will be longer than 3 characters e.g 1p21 so this is the simplest filter
#I can put into block chromosome or arm only annotations.
      next;
    }

    #Occassionally need to do some tidying up of the chromosome location
    my $cloc;
    if ( $data->[1] =~ /(.*)-.*/ ) {
      $cloc = $1;
    } else {
      $cloc = $data->[1];
    }
    $cloc =~ s/\s*//ig;    #have to trim these things
                           #Load the feature into the composite and featureLocations hashes
    if ( ( !( $data->[2] ) ) || $data->[2] =~ /NA/i ) {

      #There is no composite here so the setsize gets incremented
      $setsize++;
      $composites{ $data->[0] }       = $data->[0];    # The cloneid is its own composite id
      $featureLocations{ $data->[0] } =
        $cloc;            #This effectively using the composite id as they are the same in this case
      $bands{$cloc}++;    #increment the band/window count
    } else {

      #A composite is present
      $composites{ $data->[0] } = $data->[2];

      #Check to see if it has not been already counted
      if ( !( exists $featureLocations{ $data->[2] } ) ) {

        #The composite has not yet been added, update and increment counts
        $featureLocations{ $data->[2] } = $cloc;
        $setsize++;
        $bands{$cloc}++;
      }
    }
  }
}

#-------------------------
#####################################################################################################
#Title:  countCompositesPerAffymetrixIDRef
#Description: Builds the reference set data for an Affymetix chip directly from Biomart using R
#Author: Jon
#Date: 21.07.2005
##############################################################################################
sub countCompositesPerAffymetrixIDRef {
  my ($affy_chip_id) = @_;
  my %comps;

  #Get the affy chip design from biomart
  #Cycle through and assign the band counts and the composites
  #increment the set and band counts
  undef %comps;
}
##############################################################################################
#-------------------------
# Title: _load_properties
# Description: written to replace hardcoding of karyotype reference files and the banding colour scheme.
# Data now held in reference file ideograms.conf in $baseRefDir
# Date: 20.07.2005
# Author: Jon
#------------------------
sub _load_properties {
  my $confile = $baseRefDir . "ideogram_conf.txt";
  open( CO, $confile ) or die "Can't open $confile\n";
  my $table;
  while (<CO>) {
    my $linedata = $_;
    if ( $linedata =~ /(\S*)\s*=\s*{/i ) {
      $table = $1;
    } elsif ( $linedata =~ /(\S*)\s*=\s*(\S*)\s*/ ) {
      if ( $table eq 'ideograms' ) {
        $ideograms{$1} = $2;
      }
      if ( $table eq 'colours' ) {
        $colours{$1} = $2;
      }
    }
  }
  close(CO);
}

#-------------------------

1;
